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1. Introduction 

The automation of lattice perturbation theory (LPT) was pioneered by Luscher and Weisz in 
Ref. [1], where lists were used to represent the information required to evaluate gluonic Feynman 
diagrams numerically. Using a similar data structure, these early results were extended in Ref. [2] 
to the case of fermionic actions, enhanced in Ref. [3] to also allow for complex smearing, and 
further refined in Refs. [4] and [5] for various fermionic actions. 

Here we propose to use a more general data structure and represent the gluonic and fermionic 
actions as well as operators and matrix elements in a symbolic manner using a computer algebra 
system (CAS). 1 This approach has found wide use in the automation of continuum perturbation 
theory, where many high- loop calculations are preformed using FORM [6] . Unfortunately, FORM 
lacks features to efficiently handle the complicated vertices that arise in lattice actions. In this 
work, we introduce a new CAS optimized for LPT. We present a flexible framework on top of said 
CAS that is capable of performing perturbative calculations using a lattice regulator as well as a 
continuum regulator and to generate contractions for non-perturbative computations. 

In the following we briefly describe the layout of the new framework and present first results 
for the one-loop 0(a) -improvement of relativistic heavy quarks (RHQ) in the Columbia formu- 
lation [7] and the matching and O(o , )-improvement of heavy-light axial vector operators in the 
on-shell limit. 

2. Layout of the framework 

In this section we give a brief overview of the new framework that consists of three C++ 
libraries: libcas, a computer algebra system, libqft, a quantum field theory library that sits on top 
of libcas, and libint, a library for the numerical evaluation of algebraic expressions of lattice loop 
integrals. 

The library libcas parses algebraic expressions from text and XML input, implements a pattern 
matching algorithm that is largely compatible with FORM, and provides further methods to manip- 
ulate and simplify general algebraic expressions. One crucial feature of libcas that allows for the 
efficient use in the context of LPT is the function map. Compared to continuum perturbation theory 
the vertices of lattice actions are complicated functions of momenta. A straightforward expansion 
of all terms that appear in a typical Feynman diagram leads to an explosion of individual terms 
that is prohibitive for LPT. This issue is addressed within libcas by the function map that allows 
to replace sub-expressions (such as vertices) by references to unique functions in the function map 
in a straightforward manner. Hence the complicated momentum-dependence is represented as a 
separate function that only has to be evaluated once, and each instance of the vertex is represented 
by a mere function call with respective momenta as arguments. We make extensive use of the func- 
tion map not only to handle the complexity of lattice vertices but also, e.g., to perform a general 
common subexpression elimination within the lattice loop integrands. 

The library libqft extracts vertices from text-representations of lattice and continuum actions 
and operators, defines common lattice actions, performs Wick contractions, and computes loop 
integrals in dimensional regularization using Passarino-Veltman reduction [8]. 

'The possibility of such an alternative approach was briefly mentioned already in Ref. [1]. 
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// create action from formula string 

from ( " sum (x) *Qb Cx ) * C " 

" + sum (il , 4) tzeta ( il ) *Ngamma (il ) *aD (il , x) + amO 

11 - sum(il ,4)*r(il)*aD(il ,il ,x)/2 11 

" + sum ( il , 4) * sum (i2 ,4) *i_*cp (il , i2 ) /4 " 

" *Nsigma(il , i2)*aF(il , i2 ,x) " 

" )*Q(x) "); 



Figure 1: The definition of the RHQ action of Eq. (3.1) within libqft is shown above. The function from 
is defined within libqft and extracts the Feynman rules from the text representation of the action. Operators 
such as the covariant derivative aD or the field-strength tensor aF are defined within libqft. Special functions 
such as sum and symbols such as the complex i_ are defined within libcas. The y-algebra is represented by 
non-commuting functions N gamma and Nsigma. For a more detailed explanation of the above code, we 
refer to Ref. [9]. 

The library libint provides a convenient framework to numerically evaluate the lattice loop 
integrands with a variety of integration methods. 

All three libraries are designed to be easily extensible and will eventually be released as open 
source; a complete list of features will be given in an upcoming publication [9]. 

3. Example: Applications to heavy quark physics 

The framework found its first applications within the heavy-quark physics project of the RBC 
and UKQCD collaborations that uses RHQ in the Columbia formulation to describe the heavy 
quarks. Recent results have been presented in Refs. [10] and [11]. Relativistic heavy quarks, 
first proposed in Ref. [12] and further refined in Refs. [13] and [7], provide an effective heavy- 
quark action for large quark masses that is smoothly connected to a fully relativistic quark action 
as the quark mass becomes small compared to the lattice cutoff. The Columbia formulation [7] 
corresponds to the lattice action 




with heavy-quark fields Q. The parameters mo, £, and cp can be tuned to remove 0{ap) discretiza- 
tion errors in on-shell quantities, where ap corresponds to the spatial momentum of the heavy quark 
in lattice units. Figure 1 shows how this action is defined within libqft from a text representation. 
If we allow for a field rotation 

Q'(x)=Q(x)+d l £ YiDiQ{x) (3.2) 

i=l,2,3 

with parameter d\, we can match the quark fields Q' to continuum fields. Such a field rotation, 
however, leaves the mass spectrum of the theory invariant. Therefore the parameters mo, £, and cp 
can be tuned non-perturbatively without knowledge of d\ ; results for the tuning of bottom quarks 
were recently presented in Ref. [14]. 

In order to test the framework, we also performed a perturbative tuning of the parameters. We 
match the bilinear 

S(p)=UQ\ P )Q\q)) (3.3) 

i 
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Figure 2: Results for the perturbative tuning of parameters cp and £ on the 24 3 ensembles of Ref. [14] for 
bottom quarks. We compare perturbative results with results of non-perturbative (NP) tuning. Perturbative 
results are given at tree level (T) or at one-loop level (1). The subscript P indicates that we use the average 
value of the plaquette for meanfield improvement, the subscript U indicates that the Landau gauge-fixed 
average link value was used for meanfield improvement. The subscript L denotes expansion in the bare 
lattice coupling, the subscript M denotes expansion in the MS coupling constant at scale 1 /a, where a is the 
lattice spacing. The estimation of the errors is described in detail in Ref. [14]. 



to the continuum for on-shell momenta and determine win from the position of the pole, £ from the 
dispersion relation, and d\ from the spinor structure. The parameter cp is obtained by matching the 
three -point function 



amp 



-1 



x(j:(Q'(p)Q'(q'))) [D{ P + q) 
to the continuum in the on-shell limit, where 



— liba 
JvjU 



0(*)fv = E<Aj(*)A*(*')) 



L<<2'(<7)eV)> {^{Q'(q)A b v (k)d{-p)) 



(3.4) 



(3.5) 



is the lattice gluon propagator with color indices a, b and Lorentz indices fl, v. Figure 2 compares 
the perturbative and non-perturbative results for cp and £ on the 24 3 lattices of Ref. [14]. For more 
details of the calculation, we refer to Refs. [14] and [9]. Within estimated errors, the perturbative 
and non-perturbative results agree. 

The corresponding C++ code for the calculation of the necessary diagrams is shown in Fig. 3, 
and the resulting diagrams for the vertex are displayed in Fig. 4. The one-loop integral for the 
propagator is evaluated to 10 -3 relative accuracy on a regular desktop computer in less than ten 
seconds. 

In Ref. [11] preliminary results for fs and fp s were presented with light domain- wall fermions. 
The decay constants are extracted from matrix elements of 



a cont 
^0 



(x) 



p bl \/Z>*Z» 



q{x)m 5 Q(x) + c 1 £ q{x)YoY5YiDiQ(x) 



i=l,2,3 



(3.6) 
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Context c ; 

// use rhq + gauge action 
ActionRHQ rhq(&c , "Q" ) ; 
ActionGAUGE gauge(&c); 

// define field rotations 
c . coefficients « "dlFT"; 
const char* QimpD — 

"(1 + sum (i , 4) *dlFT (i ) *Ngamma (i) *aD (i , x) ) *Q (x) " ; 
FieldRotationRHQ Qimp(&c,"Q" , " QimpmomT " , QimpD); 
const char* QbimpD — 

"Qb(x)*(l - sum(j ,4)*dlFT(j)* Ngamma (j)*aDl(j ,x)) " ; 
FieldRotationRHQ Qbimp(&c , 11 Qb " , " QbimpmomT " , QbimpD); 

// perform wick contractions 
Wick w(&c ) ; 

w « rhq « gauge « Qimp « Qbimp ; 
Expression* vertex — w. contract ( 

" sum (X , mom) * QimpmomT Cq)*aACmom(mul,al,k)*QbimpmomT(-p)" ,3); 
Expression* prop — w. contract ( 

"sum(q,mom)* QimpmomT Cp) * QbimpmomT Cq) " . 2 ) ; 



Figure 3: The results of Fig. 2 are obtained using the code displayed above. The code defines the field 
rotations of Eq. (3.2) and generates the one-loop diagrams for the propagator of Eq. (3.3) and the vertex of 
Eq. (3.4). The resulting diagrams for the vertex are shown in Fig. 4. A more detailed explanation of the code 
shown above is given in Ref. [9]. 




Figure 4: One-loop diagrams of the quark-gluon vertex of the RHQ action in the Columbia formulation 
defined in Eq. (3.4). 



where A™" 1 is the temporal component of the heavy-light axial vector operator in the continuum, 
q is the domain-wall quark field, and Q is the relativistic heavy quark field. The lattice matrix 
element is related to the continuum version using the non-perturbatively determined light-light 
vector matching factor Zy and the perturbatively determined coefficients r\^i = Pm \J Z v b and c\ 
with heavy-heavy vector matching factor Z v b . The coefficient c\ is tuned to remove 0(ap) errors. 
In Tab. 1 we present preliminary meanfield-improved perturbative results at one loop for the coef- 
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(Tfa)(°> Ofag (d)(°) (cQW 

Plaquette 3.2615 -0.0099 0.0727 0.0257 
Landau link 3.3195 -0.0107 0.0753 0.0254 

Table 1: Preliminary perturbative results for rj/,/ = (n/,/)' - 1 + £ 2 (t/m)^ an d c i = ( c i)'°' ) + £ 2 (ci)( 1 \ where 
g is the strong coupling constant, for the 24 3 ensemble used in Ref. [11]. The results are given for meanfield 
improvement using the average plaquette or the Landau gauge-fixed average link. 



// use rhq + gauge action 
ActionDWF dwf (&c , " L 11 ) ; 
ActionRHQ rhq(&c , " Q " ) ; 
ActionGAUGE gauge(&c); 

// define field rotations (RHQ) 
c . coefficients « "dlFT"; 
const char* QbimpD — 

"u-(l/2)*Qb(x)*(l - sum(j , 4) * ( dlFT ( j ) /u) * Ngamma ( j ) * aDl ( j , x)) " ; 
FieldRotationRHQ Qbimp(&c , 11 Qb " , " QbimpmomT " ,QbimpD); 

// define field rotations (DWF 4d fields) 
const char* ID — 

"u" (1/2)* sum (1, 1 , N5)* (PL*delta(l , 1) + PR* delta (1 , N5) ) *L (1 , x) " ; 
FieldRotationDWF l(&c , "L" , "ImomT" ,1D ) ; 

// define improved vector operator 
c . commuting « "clFT"; 

FieldOperatorBilinear Vimp( c , " Vimp " , " lbmomT " , 
" Ngamma (rho )*( 1 + sum ( j , 4) * ( clFT ( j ) /u) *aD (j , x) *Ngamma ( j ) ) " , 
" QmomT" , 2) ; 

// perform wick contractions 
Wick w(&c ) ; 

w « rhq « dwf « gauge « 1 « Qbimp « Vimp; 
Expression* hl_op — 
w. contract ( 11 sum (P,mom)*lmomTCq)*Vimp(P)*QbimpmomTC-p)" ,2); 



Figure 5: The above code was used to obtain the results of Tab. 1. Field rotations are used to define the 
four-dimensional domain-wall quark. For a more detailed explanation, we refer to Ref. [9]. 

ficients 7]y and c\ for the 24 3 ensembles used in Ref. [11]. The perturbative results are based on 
the code displayed in Fig. 5. 

4. Concluding remarks 

A main virtue of the new framework presented in this paper is its versatility through the use 
of the new CAS. It provides a unified environment to perform perturbative calculations with both 
lattice as well as continuum regulators and is also able to generate contractions for non-perturbative 
computations. All expressions are stored in an algebraic manner and can be modified in a straight- 
forward way using pattern-matching techniques. 

We presented results of a first application in the context of the heavy quark physics program 
of the RBC and UKQCD collaborations. The framework is designed with higher-loop calculations 
in mind; a challenge that we plan to address soon. A publication containing further details is in 
progress [9]. 
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